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CHAPTER 1 
THE SCREENED FIELD OF A TEST PARTICLE 



Robert L. Dewar 

Research School of Physics & Engineering 
The Australian National University 

Cenberra ACT 0200, Australia 
E-mail: robert.dewar@anu.edu.au 

The screened field (forward field and wake) of a test particle moving at 
constant velocity through an unmagnetized coUisionless plasma is calcu- 
lated analytically and numerically. This paper is based on unpublished 
material from my MSc thesis, supervised by the late Dr K. C. Hines. 

1. Introduction 

Interest in the kinetic theory of interacting charged particles at University 
of Melbourne developed from the work of Dr Ken Hine^ in the 1950s — 
improving on a calculation by Landau^ and using a formalism developed by 
Fano,^ he calculated the slowing down distribution function of a charged 
"test particle" passing through a thin layer of material. The plasma theory 
group developing around Ken attracted a number of research students, 
including myself. The test particle diffusion problem provided a focus for 
a reading group on statistical physics^ and plasma kinetic theory,^ ^ and 
research on the problem developed in several ways, including a relativistic 
generalization. ^ 

The friction and diffusion coefficients of the Fokker-Planck equation 
suffer logarithmic divergences arising from the long-range nature of the 
Coulomb interaction. The fusion plasma theorists of the 1950s handled 
these divergences by the rather crude device of cutting off the interac- 
tion at the Debye length l/fcoi the characteristic length occurring in the 
screened electrostatic potential, gexp(— /cDr)/r, in the neighbourhood of a 
static particle embedded in an ionized medium. While the logarithmic na- 
ture of the divergence makes the coefficients only weakly dependent on the 
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cutofFplthe ad hoc nature of this procedure was not very satisfactory from 
a fundamental point of view and this sparked research in the international 
theoretical physics community to find better approaches. 

In 1960, using very sophisticated formalisms, BalescrP- and Lenard^ de- 
rived kinetic equations of the Landau form in which dynamical screening 
was incorporated through a frequency (w) and wavenumber (k) depen- 
dent dielectric constant e(w, k). The same Fokker-Planck-Landau equa- 
tion (which is now known as the Balescu-Lenard equation) was derived 
by Thompson and HubbarcCHinnil from simpler statistical physics argu- 
ments in which the diffusion coefficient was calculated from a fluctuation 
spectrum obtained by superimposing the dielectrically screened fields of 
independently moving particles. 

This approach was called the dressed test particle picture by Rostoker.'^^ 
In this approach the unperturbed "test particles" replace the actual parti- 
cles in the plasma. (In reality the trajectories are perturbed slightly by the 
fluctuations, giving rise to the linear dielectric response.) 

The dressed test particle picture was developed in Fourier, (a;,k), rep- 
resentation rather than in real space-time, (x, i), and thus it was difficult 
to visualize the actual nature of the screened potential surrounding each 
particle. In Chapters 2 and 3 (unpublished) of my Master's Thesis which 
was supervised by Ken Hines, I calculated the screened potential in real 
space for a nonrelativistic plasmaQ This contribution to the K. C. Hines 
memorial volume is based on those chapters with only a few changes to 
make it self-contained and to improve readability. 

For zero test particle velocity the solution is just the well-known De- 
bye potential exp(—A:D ?')/'', but as the velocity is increased we may expect 
qualitative changes to occur. The work of Pines and Bohm^^ in which they 
considered forced vibrations of the collective coordinates indicates that for 
a particle moving slowly with respect to the electron thermal velocity the 
Debye potential is distorted into a set of spheroidal equipotentials centred 
on the particle and still decays exponentially with distance. We show that 
neither result is true, since the criterion fc < fco is not an adequate criterion 
for the existence of collective coordinates. Similar results to ours were found 
in this case by Randilli^ by considering individual particle trajectories in 
the self consistent field. This is simply a way of circumventing the use of 



'^Chapter 1 contained a covariant rclativistic plasma response function formalism that 
was later incorporated into a paper on energy-momentum tensors for dispersive electro- 
magnetic waves.CS 
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the Vlasov equation but is equivalent to it. Rand makes an analogy be- 
tween the symmetrical result of Pines and Bohni and the Gibbs paradox of 
fluid dynamics, and this analogy is upheld by the fact that Majumdaill3l2l 
and CoherPS get results in agreement with Pines and Bohm by using fluid 
dynamical treatments. 

In this paper we consider a homogeneous isotropic magnetic-field-free 
plasma. We do not however restrict ourselves to a one-component plasma. 
In connection with the field around a small satellite, considered by Kraus 
and Watsonji^ it is necessary to consider the ions since the satellite velocity 
is comparable more with the ion velocity than the electron velocity. In this 
plasma there are two modes of longitudinal excitation, namely, ion acoustic 
waves and electron plasma waves, but if the test particle is much slower than 
the electron thermal velocity then only the ion waves can be excited — the 
case considered by the above authors. For ion waves not to be Landau 
damped out of existence the electron Debye length must be much greater 
than that of the ions, implying the electrons are much hotter than the 
ions or that the ions have a much greater charge. Kraus and WatsorP^l also 
consider the case of a dense plasma, in which local thermal equilibrium may 
be assumed. 

Despite the preceding remarks about two-component plasmas the case 
of infinite electron Debye length, which essentially reduces the problem to a 
one component model, has received the bulk of our attention. Most of our re- 
sults therefore can only directly be compared with those of MajumdarEHMl 
We show that his neglect of Landau damping is not justified near the edge 
of the wake. The supersonic case receives brief attention in Sec. |3.3| and |3.5| 

We note that Pappert^ considered the effect of an ambient magnetic 
field on the wake of a test particle, but we do not treat this case. A con- 
siderable body of Russian work on the details of the satellite problem had 
also appeared by the time of this thesisBll^^P^P^I 

2. The Formal Solution 

The electrostatic potential (^(x, t) in the vicinity of a test particle of charge 
q moving rectilinearly at velocity Vq (uo <C c) through a homogeneous, sta- 
ble dispersive dielectric medium is obtained from standard linear response 
theory as 

q f cPk exp[ik-(x ~ vot)] 
"^^ToJ (27r)3 Fe(k.Vo,k) 



(2.1) 
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(Henceforth we consider the time to be t = 0, when the test particle is at 
X = 0, or, equivalently, represent the potential in a frame moving with the 
test particle.) 



Im fD . 










co/k =jO 


i 





Fig. 1. The locus of <I>(cj/fc) in the complex plane as the phase speed, uj/k, is varied 
from — oo to oo. The point corresponding to Lu/k = is the intersection of the locus with 
the Re<l? axis as indicated. 



In Eq. (2.11 e(k-vo,k) is the dielectric constant, which, for an isotropic, 
coUisionless unmagnetized plasma is given by 



e(t.,k) = e(u;,|k|) = l 
The function $ is defined by 



$(w/fc) 

fc2 



/OO / / 

dv 
-OO 



uj/k — V -\- iO 



(2.2) 



(2.3) 



LOps denoting the plasma frequency, (Cgns/eo^s)^^^ (SI units) for species s, 
with the unperturbed number density, mg the mass, and gsiv) the one- 
dimensional projection of the unperturbed velocity distribution function 
/s(v). That is, in an arbitrary x,y,z Cartesian coordinate system, 



9s{Vz) 



1 



dvxdvyfs{\r) , 



(2.4) 



the normalization factor l/n^ being introduced so that / dvQg = 1. Hence- 
forth we take s to denote electrons and a single species of ion, denoted by 
subscripts e and i respectively. 

The assumed form of the "hodograph" of $ for each species is sketched 
in Fig. [T] which is such as to ensure stability towards exponentially growing 
oscillations by the Nyquist criterion [i.e. that the hodograph for e(w, k) not 
enclose the origin]. 
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Equation (2.1) becomes now 



= ^ [ c^^fc exp(^k^x) 
eoJ (27r)3 fc2 + $(k.vo) ' 

This is the starting point for all the following calculations. We shall con- 
sider only test particle velocities much less than the mean electron velocity 
so that we may approximate $c(w/fc) to the static value /c^g = $e(0), the 
square of the (generalized) inverse Debye length for the electrons. Thus 

Hu/k) « fcL + Wpi r dv ,/'^'"\ . (2.6) 



uj /k — V + iQ 

Note that if the electrons are extremely hot with respect to the ions 
then fcp^ — ^ and the ions form an essentially one component plasma with 
the electrons forming a neutralizing background. The physical explanation 
of this is that, although the electrons are much lighter than the ions, they 
are moving too fast to be appreciably deflected from their paths within the 
range of the test particle and hence do not take part in the screening. If, on 
the other hand, the test particle has a velocity comparable with that of an 
average electron then the ions may be regarded as a uniform background 
due to their inertia. This case is therefore formally identical to the case 
4o = 0. 

We give finally the normalized Maxwellian and Lorentzian distribution 
functions for a species s in a non-relativistic plasma together with the cor- 
responding polarization functions $s. 

Maxwellian case: 

^^^"^^fej ^^P-^=(^^)V^^^P-2(,^j ' ^'-'^ 
where nig is the particle mass and Tg the temperature in energy units, giving 



1/2 / 



+ '[2) "^"^Py 2 



(2.8) 



for any dimensionless x, where kf^^ = e'^Us/enTs and 

^{y) = dt expt'^ . (2.9) 



Lorentzian 



case: 







9siv) = -^^. (2.10) 
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fr2 



(2.11) 



where kus = t^ps/^so- 



3. Analytical approximations 




Fig. 2. The polarization function ^{uj/k) (in units such that Wpi = fcjj; = 1) for a 



Maxwelhan ion distribution and electron Debye constant 



0.5 (dashed horizontal 



line). The asymptotic expansion fc^^ — (fc^/w^ + v^k"^ /u}'^) is included for comparison, 
being represented by the dashed curve. 

Although it is clearly impossible analytically to evaluate the integral 
Eq. (2.5) for general g{v), or even for as simple a distribution as the 



Lorentzian, one may derive approximations for the integral in various ranges 
of Wo that enable one to gain some understanding of its properties. 
First we sketch (^{uj/k) and its asymptotic expansion in Fig. [2j 



3.1. Small Vq 



We note from equation Eq. (2.5) that the argument of $ = $o + 'I'i is less 



than or equal to Vq. Consquently for small only the behaviour of $ near 
the origin will affect the potential. 
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Expanding about the origin, <^{w) = fc^ + iaQW — aiw"^ + 0{w^), where 
oo and ai are positive constants depending on the distrubution function, 
we get, in units such that fc^ = fc^j + fc^^ = 1, 



1 1 iaow 

PT¥h ^ + 1 " + 1)2 

ai 



Then 



V(x) 



(Pk 
So J (2^' 



(fc2 + l)2 (fc2 + l)3 

k-vo 



(3.1) 



(fc2 + 1)2 



(7 I cxp(— r) 



1 

ai a-o 

(fc2 + l)2 (fc2 + l)3 

d^fc x-k 



vo-kk-vo 



47rr 



VoVo 



laovo-x 



(27r)3 (fc2 + 1)2 



(P + 1)2 (fc2 + 1)3 



(3.2) 



(27r)3 

where • denotes a unit vector. 

The first integral in the second line of Eq. (3.2) may be evaluated as 
follows 



d^k 



(27r)3 (fc2 + 1)' 



1 



(27r)2 7_/^y„ ^^^(^2 + 1)2 

d 



1 d /I 

rRc 



i(27r)2 dr \r dr 



■q{r) , 



where ri{r) is defined as follows 



dx 



exp(ix) 

x2 + z2 



(3.3) 



(3.4) 



for I argz| < 71/2 and by analytic continuation elsewhere, cutting the com- 
plex plane along the negative imaginary axis. For details of the properties 
of rj see Appendix I of the thesis.Ii^ 



The tensor term in the second line of Eq. (3.2 1 may be evaluated by 
contour integration, yielding 



q \ exp(-r) x-vq ^ d ( I d . 

+ vo- [7||xx + 71(l-xx)] -voj , 



(3.5) 
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where 



"0 



2 4 4 \ ^ 4 



^ 8 16 16 



16 



1 

8^ 





■ 2 




2 










1 










r 





IH h^ + ^|e 



(3.6) 



(3.7) 



It is to be noted that, contrary to appearances, 7^ and 7j_ are regular at 
r = 0. 

For large r, the asymptotic form of ri[r) ^ l/r holds and the exponential 
terms may be neglected. Thus, for r — >■ oo, 

q faox-vo (ag - ai) (x-Vo)2 



<^(x) 



£0 I 27r r-^ 



47r 7''^ 



(3.8) 



The first term is dominant for |x-vo| 3> i.e. for all directions not 
oblique to the direction of motion. The field is reminiscent of a dipole 
field except that it decays more rapidly with r. This may qualitatively be 
interpreted as meaning that the centre of the screening cloud has been 
displaced to a position behind the particle. It is this asymmetry that gives 
rise to the drag on a very heavy particle moving at subthermal speeds and, 
since this drag is the summation of the field of each individual particle 
moving in the average field of all the others, one supposes that the drag 
thus calculated includes single particle effects to the extent of validity of 
the linearized theory]^ i.e. except for very close collisions. 

RancP^l has calculated the screened field to first order in vq by consid- 
ering the trajectories of individual particles, explicitly rather than through 
the Vlasov equation. He obtains a function (t>i{x) as a triple integral, which 
he apparently evaluates numerically but which agrees exactljj^ with values 
calculated from ^Kc{d/dx){d/dx — l/x)rj{x). 



''This is confirmed by the work of HubbarcfHl who shows that the friction coefficient is 
the sum of the "self field" term and a fluctuating microfield term which vanishes in the 
limit of infinite test-particle mass, as do the higher transition moments. 
'^There is an error in Eq. (26) of Ref. 2, namely that the second Z in the equation should 
be replaced by 2Z/(Z + 1). 
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It is interesting to note that the sohition we have obtained differs com- 
pletely from that of Pines and Bohm^^ and workers who have used lin- 
earized fluid dynamical equations EHUI We ascribe this discrepancy to these 
authors' use of a collective approach in a region in which it is not valid: In- 
spection of Fig.|2]shows that the region in which collective coordinates exist 
can only be |w/fc| > u?, which is never entered in the case of a slow particle. 
The criterion k < /cdi is not adequate in the case of forced oscillations since 
ui may be much less than Wpi. 

We note that behind the test particle the potential has opposite sign to 
its unscreened value, which we refer to as positive. At right angles to the 
direction of motion the first order term vanishes and the second-order term 
dominates. We see that (^(x) goes negative in this direction if Og > oi. For 



the Lorentzian distribution we have from Eq. (2.111 ^i{x) = l + 2ix — 3x^ + 
^Dc + 0{x^), in units where cjpj = k^i = 1. Thus gq = 2 and ai = 3, and 
therefore a§ > oi so that the 90° potential does go negative somewhere 
between r = and r ~ oo. 

Wc shall now discuss the situation when the particle has sufficient ve- 
locity for wave excitation to be possible. 

The dispersion relation for plasma waves is e(cj, k) — 0, i.e. 

<i>{uj/k) = -k^ , (3.9) 

where uj is in general complex, but if the waves are but slightly damped 
then we may take uj real as a first approximation. It is then clear that only 
the regions in which Re<I'(a;/fc) is negative are available for excitation. In 
Fig. [2] this means the region a < uj/k < b which, it will be noted, exists if 
/cQg is sufficiently small. 

This is the region of ion acoustic waves, which resemble sound waves in 
ordinary gases in that there is an upper bound, 6, to their phase velocity. To 
guide our mathematical investigations we shall endeavour to give a physical 
picture of of the processes that can occur. 

If a particle has a velocity vq, a < vq < b, then it can excite waves with 
the same phase velocity as its own velocity. We might therefore expect 
it to be followed by a train of waves with this phase velocity, spreading 
out laterally with the maximum group velocity possible. This is roughly 
indicated in Fig. [3] 

If ti > 6 then we can expect no monochromatic train but some kind 
of shock wave confined within the "Mach cone". This is analogous to the 
formation of the shock wave behind a supersonic object in the atmosphere 
or to a longitudinal Cerenkov radiation. 
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velocity 



phase _ 
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Fig. 3. Schematic of wake field behind a charged test particle moving through a plasma. 



A note of caution must be sounded at this point with respect to the 
spreading out of the wake. Since this is controlled by the interference of 
the generated waves it is critically dependent on the dispersive nature of 
the medium, as will be discussed later. It is well known that in regions of 
anomalous dispersion in optical media the concept of group velocity breaks 
down entirely.^^ We thus expect the situation to be complicated if Landau 
damping becomes important. 

Since we are seeking wave behaviour a contour integral method is the 
obvious choice and has been much used, see e.g. Ref. [TBI We indicate the 
general method below as well as some of its analytical difficulties, not men- 
tioned in the literature due in effect to the assumption from the outset that 
an asymptotic form for the dielectric constant is a valid approximation. 



In equation Eq. (2.5) use cylindrical coordinates with axis parallel to 
vq. Then 



q f°° 
V5 = -„ / dfc_L k^Jo{k^x±_) X 
(27r)^£o Jo 



dki exp(iA:ia;i) 



kivo 



Xkl + kl) 



1/2 



(3.10) 



The contour in the ki integral may be completed in the lower/upper half 
of the complex fci plane according as Xi ^ 0. However, note that, since ki = 
±A;_L are branch points, the contours must be indented as shown in Fig. |4] 
The zeros of the denominator will give rise to poles whose contributions 
may be evaluated by the method of residues. 

For reasonably stable distributions, at least, one may show that the 



plasma wave poles corresponding to Eq. (3.9) are in the lower half fci plane. 



Consequently there is a wave excitation only behind the particle. There are 
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Fig. 4. Contours used in the evaluation of the ki integral in Eq. | |3.10| . 



other zeros of e somewhere in the lower half plane, but since they give rise 
to strongly damped contributions, one hopes that they may be ignored, at 
any rate for large |a;i|. Indeed the contribution from the region of the cut 
will also be negligible compared with the wave contribution. 

For xi > it is necessary to consider the contribution from the region 
of the cut. The part nearest the real axis is clearly the most important for 
large so let us consider the region around ikj_ 

In this region \kiVo/{k1 + kj_y^'^\ is large and the asymptotic expansion 



for $ may be used. This may be obtained from Eq. (2.6 1 by expanding the 
denominator and integrating term by term. 



fr2 



27nei~lmz)g[{z)Ljl, 



(3.11) 



where the last term is required for analytical continuation into the lower half 
z plane. Here vf is the mean-square ion speed (the square of the ion thermal 
speed in a Maxwellian plasma). However, if ki is in the lower/upper half 

1/2 

plane, so also is ki/ (fc^ + k^} ■ Hence the analytical continuation term is 
not required in the region we are considering. At ki = —ikj_ it gives rise in 
general to an essential singularity, but this contribution is being neglected 
as mentioned above. We now note the curious fact that the asymptotic 
expansion is single- valued in the upper half plane since only even powers of 

1/2 

(^ki + fc^) occur. Thus there is only a simple pole in the neighbourhood 
of fci = ik±. 

For simplicity let us retain only terms to order 1/z'^. Note that to this 
order the maximum phase velocity is 



b = Wpi/ZcDe = Cs , 



(3.12) 
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the constant Cg being commonly called the ion sound speed. 



3.2. Intermediate velocity (vf)^'''^ <^ Vq < Cg (forward field) 

If vq <ti b then the pole is approximately at fci = ik± and we may take 
^Dc = without greatly altering the situation. Then 



q 

47reo 



«o Jo 



°° Jo {uJpiXj_y/vo) exp {-ujpiXiy/vo) 

dy 



1 



y 



(3.13) 



An alternative form is obtained by noting thalP^ 

1 f°° 
= dk± Jo{k±x±) exp{~xik±) 

''' Jq 



"pi 



dyJo . 
vo Jo V 



/ 0JpiX±y \ 



Thus 



Atteo vo Jo 



dy 



1 



y 



;Jo 



Vo J 



^^^\x\_y 

Vo 



cxp 



^^iXiy 



exp 



i^pjXiy 
Vo 



(3.14) 



(3.15) 



Note that the characteristic length is now vo/ujpi, the distance the par- 
ticle moves in an ion plasma oscillation time, rather than the ion Debye 
length. The new length is longer, and furthermore, the screening is again 
not exponential at large r as we shall now show. 



For xj_ = : ip = 



1 



47reo xi 



1 vrj 



Vo 



lUlpjXi 

Vo 



q 2ujpi f LOpiT 



For xi — : ip 



Atteo Vo 

_q L 

47reo x± 



Vo 



-3 



UJpiX_l_ 



Vo 



lo 



iiJpiX±_ 
Vo 



(3.16) 



Lo 



UJpiX±^ 

Vo 



Attso Vo V ^'0 

where Iq and Lo modified Bessel and Struve functions respectively. 



(3.17) 



A general asymptotic form at large r may be obtained from Eq. (3.15) 



by approximating the denominator to 1 and differentiating Eq. (3.14) twice 
with respect to WpiXi/wo. Thus 

q 2(x-vo)2 - (xxvo)^ 



(p - 



47reo Vo (wpir/uo)^ 



(3.18) 
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Note that this obeys the same inverse third power law as in the very 
low velocity case, and also as in this case the field goes negative, for large 
r, at large angles between x and vq. The angle at which the field changes 
sign is approximately 55°. 

3.3. Supersonic velocities Vq > Cg (forward field) 

As Vq increases beyond maximum phase velocity b (the ion sound speed) 
the pole moves from ik± to i/cdc and the forward field changes to a Debye 
potential with Debye length k^^. This has the physical interpretation that 
if the particle velocity greatly exceeds the maximum ion wave phase velocity 
then the ions are too sluggish to participate in the screening. 

We must now discuss the case xi < where plasma wave excitation is 
assumed to be the dominant contribution. The positions of the plasma wave 



poles are given by the dispersion relation Eq. (3.9). It will be assumed that 



the asymptotic expansion Eq. (3.11 ) represents a valid approximation to 4>, 
with the one alteration that, since z is on or near the real axis, the last term 
is —i'Kg[{z). [It will be noted that in a sector containing the real axis g[{z) 
usually decays exponentially, and the factor multiplying it is irrelevant to 
the asymptotic expansion; but we require an approximation to $ and so 
analytically continue off the real axis where the in factor is known to be 
exact.] 

The dispersion relation is now 

k' + kl^ - + - (^) = . (3.19) 

This is valid for \Lo/k\ ^ (uj^)^/^, i.e. when the "finite temperature" 
term vfk'^/w'^ is small compared with fc^/cj^. 

3.4. Intermediate velocity (vf)^^^ < Vq < Cg (wake) 

We again make the simplifying assumption that fcoc — 0. The dispersion 
relation can now be reduced to the familiar form 

,2 ^ , ,2 I ~;2l2 I ,2 /'^Pi^ ^ 



, k' + zTTo^^iSgn^ g'l^ =0, (3.20) 



To simplify the analysis let us use units such that /coi = Wp; = 1. In 
these units vf — 3 (Maxwellian case), and generally we may assume vf ^ 1. 

The finite temperature correction is small provided fc^ 1, i.e. for 
wavelengths much greater than a Debye length. Thus, for the dispersion 
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relation to be valid, we must require both |fci| ^ 1 and fcj^ <C 1. Since 
(jj = fci^o ~ 1 it is clear that we require 3> 1. 

Substituting lu — kiVQ and — kf + k'j_, and assuming k± small, we 
find 



kl = ±- 



1 

■^0 



1 



2 I 2 

avQ kj_ 



±- 



'1 + a 






'1 + a 


^kl 


"")] 










L "^o 













(3.21) 



where a = vf /vq <C 1. 

The last term takes account of Landau damping and, being nega- 
tive, serves to displace the poles in Fig. |4]below the real fci axis. Neglecting 
the damping term, which is very small until kj_ gets close to 1, we find 



d 

dki 



fc2 + fci + $ 



kiv 

Skl + kjf'J\ 
2 



polos 



= ± 



where y = k±vo. 

The two factors containing a are small correction terms. Using 



Eq. (|3.21| to locate the poles and Eq. (|3.22|) to evaluate the residues we 

(3.23) 



find from Eq. (|3.10| that 

2-a 



f ■ 
where 



•^0 









Xi 


)'{! 


Xi 




exp ^— i 








'^'2j. 










Vo 



Im 



2 (1 + a + y2) ^ [^{l + a + y2)i/2 ) ^^''^^^ 

is the Landau damping term, which has the effect of rapidly cutting off the 
integral when \y\ gets close to wq, and / is the integral 

Joiby) cxp(-ia?/^'' 



Iia,b,c) = I dyy 

lo (l-|-r)(l + cy2) 



ith 



^(y) = exp 



7(2/) 



(3.25) 
(3.26) 



In Eq. (3.231 the subsitutions for the dummy arguments a, b and c are 

^ _ 

Vo ' Vo 



a 

"=2 



x± a 
— and c — — 



(3.27) 
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When 5 = we suppose that the cutoff factor ?/' may be neglected. Then 
/(a,0,c)«^^^ e^^E,ita)~e^^/^E,(^-^^ (3.28) 

where Ei is the exponential integral.!^ 

To obtain an asymptotic form in the general case it is helpful to trans- 
form the one-sided integral into a two-sided integral by a method due to 
Hankelj^ noting that, for a; > 0, 

1 



(3.29) 



where Hq^\z) is a Hankel function, defined on the complex z-plane cut 
along the negative real axis, being regular elsewhere. Then 



I{a,b,c) = 



1 



(l + r)(l + c?y^) 



(3.30) 



If a were zero and ip well behaved, we could evaluate the integral by 
completing the contour in the upper half plane. Instead we seek the saddle 
points and poles of the integrand, the poles being at y = ±z and y = 
±i/^/c. To find the saddle points let us suppose that the asymptotic form 
fori7W(z), 



2 expi(z — j) 



|z|^0. 



(3.31) 



is valid in the region of these points. (The asymptotic form will be valid in 
the region of the poles provided |6| 1, c < 1.) 

Since the denominator is slowly varying away from its zeros the saddle 
points are determined by the exponential factor, the exponent of which, 
being quadratic, has only one stationary point. This is at y — yo, where 
2/0 = b/2a. For Eq. (3.30) to be vahd we require \yo\ 1. 



We conclude that, provided the two conditions 6 >• 1 and 6 2|a| are 
satisfied, we may approximate / by 



I{a,b,c) w ^— •\/27r 



Vb 



dy jy + iO)^/^ e^p{~iay'^)tlj{y) ^^^^ 

(l + y2)(i + cy2) • ^^-^^^ 



This is a Fourier transform, which may be transformed by the convolu- 
tion theorem to 



/(a,6,c)«^^^,* 



° {y + ^0)'/'exp{~^ay^)e^''y 

"■y rr, — — 2\ ' 3.33) 



August 2, 2011 5:1 



WSPC/Trim Size: 9in x 6in for Review Volume 



HinesMem' Dewar 



16 



R.L. Dewar 



where tpb = 'ip{y)e'^^^dy /2it is a delta-like function of b with width I/vq, 
which will have the effect of destroying all wave structure with wavelength 
of a Debye length or less. 

The integral may now be estimated by deforming the contour of integra- 
tion to cross the saddle point along the line of steepest descent. (The pole 
at ij \fci being far up the imaginary axis, gives negligible contribution.) 




Re y 



Fig. 5. The contour used for estimating the integral in Eq. l |3.33| l. 



We make the further requirement that the width of the saddle at ?/o be 
much less than yo- That is, l/-\/a ^ j/o- Then 



/(a, 6, c) w -06 * 



yy^exp(7;fe2/4a) 




(3.34) 



Since the first term varies slowly compared with i/'fc, the convolution will 
leave it but little changed, whereas the rapid fluctuations in exp(i&^/4a) 
will be damped out. Supposing that the exponential term may locally be 
approximated by a monochromatic wave, we finally obtain 

T e''-^ i j/q/^ exp(i&74a) 
2b 1-c 



/(a, 6, 



(3.35) 



c V2a6 (1 + %')(! + C2/2) ' 
with yo = b/2a. 

The method of steepest descent we have used to approximate the inte- 
gral is equivalent, in the neighbourhood of the saddle point, to the method 
of stationary phase employed by MajumdaJ^ but takes better account of 
the behaviour away from the saddle point by including the contribution of 
the pole aX y — i. We have also attempted to take some account of Landau 
damping, which the above author was unable to do owing to his formulation 
in terms of fluid dynamics. 
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Below we recapitulate the criteria for the validity of Eq. (3.35) 
c<l, 6>1, 6>2a, 6>2ya. 



(3.36) 



With a, b and c given by Eq. (3.271 these correspond physically (recalling 
that in this subsection we are using units such that fcoi ~ ^pi = 1) to the 
inequalities 



« 1, 



Vo 



> 1, a;j_>a|xi|, a;_L > (wf jxiD^/^ . 



(3.37) 




Fig. 6. Schematic of the wake region of an intermediate-velocity test particle indicating 
where the approximations used to derive Eq. ( |3.35| l break down (hatched region) and 
the "thermal Mach cone" (dashed line). 



Observe that the Landau damping term cuts off the second term at 
about x± ~ \xi\/vo, which might be called a thermal Mach cone. This be- 
haviour can be understood from the fact that plasma waves with high group 
velocity have low phase velocity, and hence are strongly damped. The first 
term, however, is unaffected by Landau damping and exhibits no thermal 
Mach cone, although it decays exponentially with x± . We indicate the ther- 
mal Mach cone and the region excluded by the inequalities in Eq. (3.371 in 

Fig.m 



3.5. Supersonic velocities Vq > Cg (wake field) 

This is the case considered by Kraus and WatsorP^l but their treatment 
neglects finite ion temperature effects on the dispersion relation and is con- 
sequently inadequate for finding variations over distances less than fcj^^ . We 
shall not go into detail on this case but shall indicate the approximations 
that may be made and the integral obtained. 
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The simplifying assumption made by the above authors can be repre- 
sented as 



(3.38) 



where b is here the upper bound to the phase velocity Cg by Eq. (3.12)], 



not the dummy argument used in the previous subsection. Taking units 



such that uj-o 



kDe 



1 we have 



vf < 1 



1 , f » 1 and » 1 . 



(3.39) 



Using the same methods as before, but neglecting Landau damping, we 
find the approximate form for xi < 



2 



dk I k 



± Kj_ 



(i + fci)3/' 



Jo(a:_LA;_L) sin 




(1 



/c2)l/2 



(3.40) 



Kraus and Watson in effect assume that the critical part of the integral 
for determining the large |xi|/fo behaviour of {x± small) is near k± ^ 
so that terms such as 1 + fc^ may be approximated to 1. However, it must 
be observed that there is a range of k± (between 1 and 1/vf) in which 
the argument of the sine function is but slowly varying. This suggests that, 
superimposed on the endpoint contribution, there may be a sinusoidal term. 

When x± is not small we may transform the integral into a two-sided 



integral as in Sec. 3.4 There will be a saddle point which disappears into 
the origin at the Mach cone xj_ — \xi\/vq, but the contribution from the sin- 
gularities of the integrand will be complicated by the essential singularities 
at fc_L = ±1. 



4. Numerical Solution 

The explorations of the previous section showed that some insight into the 
nature of the fields around a particle in a coUisionless plasma, as given 



formally by Eq. (2.1), may be obtained by asymptotic analyses in various 
limits. However, these are often difficult and the errors difficult to quantify. 
It is thus essential to supplement such analytic work with numerical and 
graphical methods (and vice versa). 

In this section we summarize how this was done computationally in 
my MSc thesis work^^ using the newly arrived IBM 7044 (which had just 
replaced the Australian- built computer CSIRACpS). 



August 2, 2011 5:1 



WSPC/Trim Size: 9in x 6in for Review Volume 



HinesMem' Dewar 



The Screened Field of a Test Particle 



19 



4.1. Formulation of the numerical method 



To evaluate the integral Eq. (2.1 1 numerically we must choose axes such 



that as much of the integration as possible may be done analytically while 
the remaining integrations are as tractable as possible. It is found that, 
although one may come tantalizingly close to reducing the triple integral 
to a single integral, it is in general necessary to evaluate a double integral. 




Fig. 7. Coordinate system used for numerical evaluation of the triple integral in 
Eq. 



Three systems of coordinates suggest themselves: the cylindrical coor- 
dinates of Sec. [Sj spherical polars with axis along Vg, and spherical polars 
with axis along x. The first method involves infinite integrals with inte- 
grands having a large peak near the plasma wave pole and we consequently 
reject it. The last choice appears to have the advantage over the remaining 
one that it is easier to understand the behaviour of the integrands (at least 
at low vq) and to isolate the singularity caused by the behaviour of ri{x) at 
X = 0. These axes are sketched in Fig.[7| In this coordinate system we have 

k'Vo = iivii + \J\ — [J? cos uj^ , (4-1) 



where v\\ = vq-x = uqcos?/', uj_ — |voXx| = uosinVj, and [i, = k-x 

k'^ cxp{ik^r) 



In these coordinates Eq. (|2.1|) becomes 

q 



LP 



(2^)3£0 



1 /'27r foo 

dfi / dip Re / dk 

1 Jo Jo k^ + ^{fiv\\ + ^/l — fi"^ cos (j>V±) 



(4.2) 



where r = |x| and $ is defined in Eq. (2.3). The units used were such that 
/cDi — LUpi — 1 and q — Atteq. 
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By using the fact that — <&(— x) for real x, the range of the fj, and 

integrations was reduced by half. It was found usefuF^ to define the new 



special function r](z), Eq. (3.4). Analytical properties of r](z) are discussed 
in Appendix I of the thesia^and a listing of a Fortran IV subroutine for 
its efficient evaluation is given. This was used in the code developed to 
calculate the results presented below. (More discussion of the numerical 
method is given in the thesis .El) 



In terms of 77, Eq. (|4.2|) reduces to a double integration 

'^r]{xV^ 



dx Im 



where 



1/2 



1/2 



cos (j) v±_ 



(4.3) 



(4.4) 



4.2. Small Vq results 




Fig. 8. Solid curves are the computed deviation from the Debye potential, as defined 
in the text, for a particle moving at velocity dq = 0.2a;pi/A;Di through a plasma with 
Lorentzian ion distribution and infinite-temperature electrons. Dashed curve is an ana- 
lytic approximation to the forward field from Eq. l |3.5| . 



The case where the test particle is travelling slower than the ion thermal 
speed was studied using the Lorentzian distribution function Eq. ( 2.10 1 with 
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^Dc = 0. The treatment of Sec. |3.1| indicates that the choice of distribution 
function and electron Dcbye length simply alters the relative magnitude of 
the constants fc^, oq and ai, so we expect there is little to be gained by 
study of a large variety of cases. 




Fig. 9. Equipotentials in a plasma the same as assumed for Fig. |8] but with vq = 
0.4ajpi/A;Di (the direction of vq is indicated by an arrow). The sohd equipotentials indicate 
positive (+ve) ip and the dashed lines negative (— ve). The scale in Debye lengths is 
indicated top right. 



In Fig. [8] we plot the deviation from the Debye potential along lines 
making angles tp of 0°, 60° and 90° to the direction of travel, the deviation 
being defined as [ip{x.)—ipj:,{r)]/(pQ{r), where vdI?") = (/exp(— A;Di?')/47r£o'' is 
the Debye potential and (po{r) = q/Aireor is the bare, unscreened potential. 
The approximate result at 0° obtained from the term of Eq. (3.5 1 first order 
in Vq is also plotted (dashed curve), the agreement being seen to be quite 
reasonable. 

The qualitative form of Fig. [9] is also in agreement with that predicted 
from the analytical work. Note the potential well behind the particle, be- 
yond which the field becomes attractive to charges of the same sign as the 
test particle. 
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4.3. Intermediate Vq results 




Fig. 10. Plots of ip{r, ip) /(po{r) as a function of distance r (in units of i;o/t*^pi) calculated 
for a test particle moving at vq = Swpi/fcoi in a plasma with Lorentzian ions and fcoo = 0. 
The forward field is labelled 0° and the wake field is labelled 180°. The dashed line is 
the unscreened case, ip = ipo- 



The intermediate case, where the test particle moves faster than the ion 
thermal speed but less than the ion sound speed Cs, was studied in part us- 
ing both Lorentzian and Maxwellian distribution functions as these to some 
extent represent extremes of Landau damping. The Lorentzian distribution 
is sufficiently simple for an exact dispersion relation to be derivable, thus 
obviating the use of asymptotic expansions of $. Again we took fcDo = 0, 
so = 00 by Eq. ( |3l2 |. 

Potentials on the axis of symmetry, parallel (forward field, ip — 0) or an- 
tiparallel (wake field, -0 = tt) to the direction of motion of the test particle, 
are plotted in Figs. [ToHT2| 



Forward field: In all three of the plots the results agreed well with Eq. (3.16 1, 



with agreement particularly good for the Maxwellian case (four significant 
figures at vq = Sajpi/fcoj. [The potentials at 90° to the direction of motion 



also agreed reasonably with Eq. (3.17) 



Wake field, Lorentzian case, Fig. \10^ Within the wake Landau damping 
had a dominant effect in the Lorentzian case even at the centre of the 
wake. This is due to the slow decay of the tail of the distribution function. 
The analytical prediction (which was not derived in Sec. [3]) was virtually 
indistinguishable on the scale of the graph, which is to be expected because 
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the dispersion relation is exact. 

Wake field, Maxwellian case: For a Maxwellian distribution function the 



agreement with Eq. (3.28) was quite good in the Vq — 8 case, Fig. 11 



but not so good in the Vq — 4 case, Fig. [12] This may be explained by the 
inadequacy of the asymptotic expansion Eq. (3.11 ) when \io/k\ < 4. Landau 
damping may also have some effect. 

Owing to computing time limitations the complete wake structure was 
mapped only in the vq — A Maxwellian case. A contour plot is given in 



August 2, 2011 5:1 



WSPC/Trim Size: 9in x 6in for Review Volume 



HinesMem' Dewar 



24 R.L. Dewar 




Fig. 13. Contours of (f{x.)/ipo{r) for a test particle moving with hq = Aujpi/kj^i (the 
direction of vq is indicated by an arrow) in a plasma with Maxwellian ions and fcoo = 0. 
The scale in units of vg/uipi is indicated lower right, the "thermal Mach cone" x±_ = 
\xi luipi/kuivo is indicated by dash-dotted lines, and the boundary of the region of validity 
of Eq. | |3.35[ l is indicated by the dotted curve (cf. Fig.|6]l. 

Fig. [13) which represented almost two hours of computing time on the IBM 
7044 and is still not very accurate. 



4.4. Supersonic case 

The supersonic case, vq > Cs, considered in Sec. |3.5| was only briefly studied 
numerically. Figure [14] shows a fairly sharp cutoff near the Mach cone but 
no shock front. Indeed the function decays monotonically off the axis in 
complete contrast with the Kraus-Watson^^ result, which predicts an initial 
increase. 
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Fig. 14. Plot of ip(r, 4>)/'Po(r) vs. ip (with r = 48/fc]3i) in the wake of a supersonic test 
particle moving at vq = 32i.jpi/fcj3i = 4Cs in a plasma with kj^^ = k^i/S. Surprisingly, 
no shock wave is seen at the expected position of the Mach cone, which is calculated 
from a;x/|a::i| = i^^pi/vokj^^ and indicated by the arrow marked M at i/" = 166°. 



4.5. Conclusion 

In the intervening years since my MSc work other authors have made similar 
calculations. I have not attempted a complete literature search, but note 
that the research project must have been topical at the time as the inverse 



third power asymptotic behaviour in Eq. (3.8) was announced a year later 
by Montgomery, Joyce and Sugihara,!^ in a paper that has been cited 58 
times. In fact the problem is even more topical today with the rise of interest 
in dusty plasmas — for instance the paper by Ishihara and VladimiroA^^Sl on 
the wake potential of a dust grain in a plasma with ion flow has been cited 
more than 80 times. 
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